Rheology of a supercooled polymer melt near an oscillating plate: an application of 

multiscale modeling 
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The behavior of a supercooled polymer melt composed of short chains with ten beads near an 
oscillating plate are simulated by using a hybrid simulation of molecular dynamics (MD) and com- 
putational fluid dynamics (CFD). In the method, the macroscopic dynamics are solved by using 
CFD, but, instead of using any constitutive equations, a local stress is calculated by using a non- 
equilibrium MD simulation associated at each lattice node in the CFD calculation. It is seen that 
the local rheology of the melt varies considerably in a thin viscous diffusion layer that arises near 
an oscillating plate. It is also found that the local rheology of the melt is divided into the three 
different regimes, i.e., the viscous fluid, viscoelastic liquid, and viscoelastic solid regimes, according 
to the local Deborah number De, which is defined with the Rouse or a relaxation time, tr or r Q , 
and the angular frequency of the plate u as De R =u>TR or De a =0JT a - The melt behaves as a viscous 
fluid when De R < 1, and the crossover between the liquid-like and solid- like regime takes place 
around De a ~ 1. 

PACS numbers: 31.15.xv 46.15.-x 66.70.Hk 
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Due to the emergence of nano-science (or nano- 
technology) , a great innovation of the functional materi- 
als of soft matters is expected recently. An important 
property of soft matter is that there exit strong cor- 
relations between the nano-scale structure and kinetics 
and the macroscopic behavior of the system. In other 
words, the nano-scale structures and kinetics may affect 
the macroscopic physical properties, and the macroscopic 
behavior may also affect the microscopic physics. This 
property of the coupling of the microscopic and macro- 
scopic physics is very useful for functional materials, but 
it also makes the theoretical and experimental investi- 
gations of soft matter difficult. The simulations may be 
useful for predicting the complex behaviors of soft mat- 
ter. The difficulty in the simulation of soft matter is 
the coexistence of the macroscopic and microscopic scales 
involved. The macroscopic dynamics are usually solved 
by using the continuum equation, e.g., the Navier-Stokes 
equation, with a constitutive relation. For soft matter, 
however, the constitutive relation is so complicated that 
no explicit formulas are obtained in general. In contrast, 
the molecular dynamics simulations are performed with- 
out using any constitutive relations. However, the time 
and space to be simulated are restricted to very micro- 
scopic scales. In order to overcome this difficulty, we have 
recently developed a hybrid simulation of the molecular 
dynamics (MD) and the computational fluid dynamics 
(CFD) based on a local sampling strategy, in which the 
macroscopic dynamics are solved using a CFD scheme 
but, instead of using any constitutive equations, a lo- 
cal stress is calculated by using a non-equilibrium MD 
simulation associated with each lattice node of the CFD 
computation. [l| The basic idea of the hybrid simulation 
method was put forward earlier by Kevrekidis et aL[2| 
and also by Ren and E.Q 

In previous papers, the validity of the hybrid methods 



and their efficiencies are examined intensively for viscous 
fluids without memory effects. De et al. have recently 
proposed a new hybrid method, which is called the scale 
bridging method in their paper, that can correctly repro- 
duce the memory effect of the polymeric liquid, and per- 
formed a simulation of a non-linear viscoelastic polymeric 
liquid between two parallel oscillating plates. Q They 
have also compared the results obtained by the scale 
bridging method with those obtained by a full MD simu- 
lation, thereby demonstrating validity of the method. In 
the present letter, we also model the behavior of a poly- 
mer melt between two parallel oscillating plates by using 
the hybrid method, but we focus on the complex rheol- 
ogy of a supercooled polymer melt in the viscous diffu- 
sion layer that arises near the fast oscillating plate. The 
boundary layer arises if the width between the plates is 
much larger than the thickness of viscous diffusion layer. 
Note that, in Ref. 0, the thicknesses of viscous diffusion 
layers are estimated to be comparable to the widths of 
the plates thus the boundary layers are not clearly 
seen. 

In the present problem, the macroscopic quantities are 
quite non-uniform, and two different characteristic length 
scales appear that must be resolved; one is that of a poly- 
mer chain and the other is that of a boundary layer aris- 
ing near the oscillating plate. This problem constitutes 
an important application of multiscale modeling since it 
is quite difficult to solve this problem by using a full MD 
simulation because the length of the boundary layer is 
much larger than that of a polymer chain. In the follow- 
ing, we briefly state the problem and outline the hybrid 
simulation method, and then discuss the numerical re- 
sults. Finally we give a summary of our conclusions. 

We consider a polymer melt with a uniform density po 
and a temperature Tq between two parallel plates (see 
Fig. HJa)). The upper- and lower-plate start to oscil- 
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FIG. 1: Schematics for the problem and mesh system. 



late in opposite, parallel directions with a constant fre- 
quency u>o at t=0. The polymer melt is composed of 
short chains with ten beads. All of the bead particles in- 
teract with a truncated Lennard- Jones potential defined 
by@ U LJ (r) = 4e [(cr/r) 12 - (cr/r) 6 ] + e for r < 2^ e a, 
and for r > 2 1 / 6 o\ By using the repulsive part of 
the Lennard- Jones potential only, we may prevent spa- 
tial overlap of the particles. Consecutive beads on each 
chain are connected by an anharmonic spring potential, 
U F (r) = -±k c R 2 ln [l - (r/R ) 2 ], with k c =30e/a 2 and 
i?o=l-5cr. The number density of the bead particles is 
fixed at p /m— 1 /a 3 , where m is the mass of the bead par- 
ticle. With this number density the configuration of bead 
particles becomes severely jammed at a low temperature, 
resulting in the complicated non-Newtonian viscosity and 
long-time relaxation phenomena typically seen in glassy 
polymers. In the present letter, we fix the temperature 
at T—0.2e/k, where k is the Boltzmann constant. 0] 

We assume that the macroscopic quantities are uni- 
form in the x- and z-directions, i.e., d/dx=d/dz=0. 
Then the macroscopic velocity v a is described by the 
equations, podv x /dt = da xy /dy and v y =v z =0, where 
<j a j3 is the stress tensor. Here and afterwards, the 
subscripts a, 0, and 7 represent the index in Carte- 
sian coordinates, i.e., {a,/3 ,j}={x ,y ,z} . We also assume 
the non-slip boundary condition at the oscillating plate, 
v x —vq coswrji at y=0 (where vq is a constant amplitude 
of the oscillation velocity), and the symmetric condition 
at y=H (i.e., v x (y=H+Ax/2)=-v x (y=H-Ax/2)). If the 
frequency lvq is large enough, the thin viscous boundary 
layer forms over the oscillating plate. The thickness of 
the layer is estimated, for fluids with a constant viscosity 
2/0 j as l v ~ 1 \J v^/loq. Note that the thickness of the vis- 
cous layer l v is much smaller than the width between the 
plates, l v -C H, but is usually much larger than the scale 
accessible to a full MD simulation, in which the charac- 
teristic length scale is the length of the polymer chain l p , 
l v > J P . 

The constitutive relation of the stress tensor is quite 
complicated || [9], Ho| : the temporal value of the stress 
tensor of a fluid clement depends on the previous values 
of the velocity gradients of the fluid element. The relation 
may be written in a functional form as, 

<Jap(t,x a ) = F a p[n a p(t',x' a (t'))], with t'<t, (1) 

where n a p is the velocity gradient, K a p = dv a /dxp, and 
x' a (t') represents the path line along which a fluid ele- 



ment has been moving. In the one-dimensional problem, 
however, we don't need to consider the path line of the 
convective fluid element since the macroscopic velocity is 
restricted to be only in the x-direction where the macro- 
scopic quantities are assumed to be uniform. Thus the 
stress tensor for the present problem may be written in 
a functional involving the local strain rate, 

a af3 (t,y)=F a0 [j(t' > y)}, with t' < t, (2) 

where 7 is the strain rate, -j=dv x /dy. Note that, al- 
though Eq. ^ becomes much simpler than Eq. (fT]), the 
temporal value of the local stress still depends on the 
previous values of the local strain rate. Its dependence 
is quite complicated, especially for glassy materials (for 
which explicit formulas are unknown in general). In our 
hybrid simulation, instead of using any explicit formulas 
for the constitutive relation, the local stress is generated 
by the non-equilibrium MD simulation associated with 
each local point. See Fig. [T] (b). 

We briefly explain our hybrid simulation method. We 
improve the previous hybrid method [l[ with regard to the 
time-integration scheme of the CFD and the treatment 
of the sampling duration of the MD so that the memory 
effects are correctly reproduced. The basic idea of the 
present method is the same as that of the scale bridging 
method proposed in Ref. 0. The computational domain 
[0,H+Ax/2] is divided into thirty- two slits with a con- 
stant width Ace. We use a usual finite volume method 
with a staggered arrangement, where the velocity is com- 
puted at the mesh node and the stress is computed at the 
middle of each slit. One hundred chains are confined in 
each cubic MD cell with a side length Zmd = 10ct. As for 
the time-integration scheme, we use the simple explicit 
Euler method with a small time-step size At. The lo- 
cal stress at each time step of the CFD is calculated by 
performing a non-equilibrium MD simulation according 
to the local strain rate in each MD cell. The techniques 
of the non-equilibrium MD simulation are the same as 
those in the previous paper we use the Lees-Edwards 
sheared periodic boundary condition and a Gaussian iso- 
kinetic thermostat to keep a constant temperature. In 
the present problem, however, we cannot assume a local 
equilibrium state at each time step of the CFD simulation 
since the relaxation time of the stress may become much 
longer than the time-step size of the CFD simulation (in 
which the macroscopic motions of the system should be 
resolved). In the present simulations, the simple time- 
average of the temporal stresses of the MD (averaged 
over the duration of a time-step of the CFD simulation) 
are used as the stresses at each time step of the CFD cal- 
culation without ignoring the transient time necessary for 
the MD system to be in steady state. Thus, the ratio of 
the time-step size of the CFD calculation At to the sam- 
pling duration of the MD simulation tun is A£/£md = 1- 
The final configuration of molecules obtained at each MD 
cell is memorized as the initial configuration for the MD 
cell at the next time step of the CFD. Thus we trace all 
of temporal evolutions of the microscopic configurations 



3 




Vo COS LOot 

(a) Polymer melt (b) Newtonian 



FIG. 2: The velocity profiles at w o t=24O7r+0, 0/tt=O, 1/2, 
3/4, 1, 3/2, and 7/4, for wo=2tt/256. (a) The result for the 
polymer melt and (b) that for a Newtonian fluid. The ver- 
tical axis represents the height y and the horizontal axis the 
velocity v x . 



with a microscopic time step so that the memory effects 
can be reproduced correctly. Note that, compared with a 
full MD simulation, we can achieve an efficient computa- 
tion with regard to the spatial integration by using MD 
cells that are smaller than the slit size used in the CFD 
simulation. The efficiency of the performance of our hy- 
brid simulation is represented by the ratio of the slit size 
used in the CFD model Ax to the cell size of the MD 
simulation Zmd, Ax/Imd- Hereafter, the quantities nor- 
malized by the units of length a and time T^=\Jm<j 1 /e 
are denoted with a superscript "*" . In the following sim- 
ulations, we fix the time-step size of the CFD simulation 
At, sampling duration of the MD simulation tyiD and 
time-step size of the MD simulation Ar as At*=ijj^ D =l 
and Ar* =0.001, respectively. Thus, one thousand MD 
steps are performed in each MD cell at each time step of 
the CFD computation. 

We perform the hybrid simulations for two cases: Case 
A, in which wg=27r/256 and £P=787.5, and Case B, in 
which ^=27r/1024 and iJ* = 1575. The amplitude of the 
oscillation velocity is fixed as i>q=10 in both cases. The 
widths of the slits are Ax* =25 for Case A and Ax* =50 
for Case B, and the ratios of the mesh size of the CFD 
Ax to the cell size of the MD Zmd are Ax/^md=2.5 for 
Case A and 5 for Case B. Figure [5] shows the instan- 
taneous velocity profiles over a period for (a) Case A 
and (b) the Newtonian fluid with a constant viscosity 
i/q=53, which corresponds to the dynamical viscosity of 
the model polymer melt for wo=27r/256. The dynami- 
cal viscosity is calculated via G\ju) (where Gi is the loss 
modulus defined below). It is seen that the boundary 
layer of the melt is much thinner than that of the New- 
tonian fluid. This is caused by the shear thinning; the 
local loss modulus near the boundary is, as we see be- 
low, much smaller than that for the linear regime, thus 
the thickness of the viscous diffusion layer becomes thin- 
ner than that for the Newtonian fluid. We also mea- 
sure the local viscoelastic properties in terms of the lo- 
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FIG. 3: The spatial variations of the local moduli Gi and 
G2 (lower axis) and the amplitude of the local strain 70 (up- 
per axis) for Case A. The dashed and dash-dotted lines show 
the storage modulus and loss modulus for the linear regime, 
respectively. The long-dashed line represents 70=2%. The 
left arrows on the right-side vertical axis show the positions 
where the local Deborah numbers, shown in Fig. 5, are equal 
to unity. 

cal storage modulus G\(y) and loss modulus C?2(y)- It 
should be noted that the local macroscopic quantities 
oscillate with a different phase retardation at each dif- 
ferent point. The local moduli are calculated in the fol- 
lowing way: The discrete Fourier transforms of the tem- 
poral evolutions of the strain 7, j(t,y) = Ljit' \y)dt' ', 
and shear stress a xy during the steady oscillation states 
are performed. The discrete Fourier transformations are 
written as g{ = J2n=o 9n exp[i2irnk/(N + 1)), with 
g l n =g{nAt, lAx) (n—0,...,N and i=0,...,32) , where g rep- 
resents the strain or shear stress (e.g., <?=7 or a xy ). By 
using the Fourier coefficients for the mode of the oscillat- 
ing plate fcoj ko=l+(uJo/2Tr)N, the time evolution of the 
local strain at y=y l can be expressed as a cosine function, 

7 '(i) =j l cos(ujot + 6 l ), (3) 

with 7 < = sj (7'1 )2 + (7"' fe0 )2 and S l = 

tan _1 (7" fe0 /7' fc0 ). Hereafter the superscripts "/" 
and "//" indicate the real and imaginary part of the 
discrete Fourier coefficients, respectively. The time 
evolution of the local shear stress can also be expressed 
as 

<J l xy {t) = a\ cos(w i + S l )-a 2 sm(uj Q t + S l ), (4) 

with a{ = a' ko cosS 1 + a" ko sirid 1 and a\ — a" ko cosd 1 — 
' 1 1 

a' ko sin 6 . Thus, the local storage modulus G\ and loss 
modulus G2 are written, respectively, as G\{y l )=a l l /^\ ) 
and G 2 {y l )=a\hl 

Figure [3] shows the spatial variations of the local stor- 
age modulus and loss modulus and the amplitude of the 
local strain for Case A. The shear thinning is seen near 
the oscillating plate; the strain 70 (or the strain rate 
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FIG. 4: The spatial variations of the local moduli G\ and Gi 
(lower axis) and the amplitude of the local strain 70 (upper 
axis) for Case B. See also the caption for Fig. 3. 
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FIG. 5: The local Deborah numbers (defined via the Rouse 
relaxation time tr and the a relaxation time r a ) De R = loqtr 
and De a = ur°. (a) Case A. (b) Case B. The right arrows 
on the t/-axis show the position where the Deborah numbers 
are equal to unity. 



70=^0 7o) is quite large near the oscillating plate and 
both moduli are much smaller than those for the linear 
regimes. In the close vicinity of the oscillating plate, the 
storage modulus G\ is much smaller than the loss mod- 
ulus G2, G\ <C G2. Hence, the melt behaves as a vis- 
cous fluid. The storage modulus rapidly grows with the 
distance from the oscillating plate, and the viscoelastic 
crossover occurs at y* ~ 200. Both moduli attain their 
linear values, which are shown as dashed and dot-dashed 
lines in Fig. [3l for distance that is far from the oscil- 
lating plate where the local strains are less than about 
two percent. The overall features are also consistent with 



Case B, although the crossover is not as clear as that in 
Case A (See Fig. 0}. Thus, the local rheology of the melt 
can be divided into three regimes, i.e., the viscous fluid, 
viscoelastic liquid, and viscoelastic solid regimes. These 
regimes may be also characterized by the two "local" 
Deborah numbers. One is defined by the local Rouse re- 
laxation time tr of the melt and the angular frequency of 
the plate ujq, De R =u)QTR, and the other is defined by the 
local a relaxation time r Q and the angular frequency ojq, 
De a =u! T a . Note that the local Rouse and a relaxation 
times vary according to the local strain rate 7, r = r(j). 
Figure [5] shows the spatial variation of the local Deborah 
number De R and De a , where the local relaxation times 
tji and r a are calculated by substituting the values of 
7q, which are obtained via Eq. ([3]), into the fitting func- 
tions for the relaxation times for the simple shear flows 
obtained in Ref. 7;. In Figs. [3] and |H the positions at 
which the local Deborah numbers become equal to unity 
are shown by the left arrows. It is seen that the melt be- 
haves as a viscous fluid, G2 ^> G\, for De R < 1, while the 
viscoelastic properties become pronounced for De R > 1. 
This is consistent with the rheology diagram for a model 
polymer melt obtained in Ref. It is also seen that 
the crossovers of the liquid- like regime, G2 > Gi, and the 
solid- like regime, G\ > G2, take place at De a ~ 1. 

In summary, we have numerically analyzed the behav- 
iors of a supercooled polymer melt near an oscillating 
plate by using the hybrid simulation of MD and CFD 
based on the local sampling strategy. In our simulation, 
the memories of the molecular configurations of the fluid 
elements are correctly traced at the microscopic level. 
The efficiencies of our hybrid simulations are represented 
by the ratios of the mesh size of the CFD simulation, Ax 
to the cell size of the MD simulation, ^md- These ratios 
are Ax/?md=2.5 for Case A and 5 for Case B. It is found 
that the local rheology of the melt varies considerably in 
a viscous boundary layer near an oscillating plate, so that 
three different regimes, i.e., the viscous fluid, viscoelastic 
liquid, and viscoelastic sold regimes, form over the oscil- 
lating plate. It is also found, in the viscous fluid regime, 
that the local Deborah number defined via the Rouse re- 
laxation time and the angular frequency of the plate is 
less than about unity, De R < 1. The crossover between 
the liquid-like and solid-like regimes takes place around 
the position where the local Deborah number defined via 
the a relaxation time and the angular frequency is equal 
to unity, De a ~ 1. 



[1] S. Yasuda and R. Yamamoto, Phys. Fluids 20, 113101 
(2008). 

[2] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. 

Kevrekidis, O. Runborg, and C. Theodoropoulos, Comm. 

Math. Sci. 1, 715 (2003). 
[3] W. Ren and W. E, J. Compt. Phys. 204, 1 (2005). 
[4] S. De, J. Fish, M. S. Shephard, P. Keblinski, and S. K. 



Kumar, Phys. Rev. E 74, 030801 (2006). 
[5] S. Sen, S. K. Kumar, and P. Keblinski, Macromolecules 

38, 650 (2005). 
[6] K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057 

(1990). 

[7] R. Yamamoto and A. Onuki, J. Chem. Phys. 117, 2359 
(2002). 



■5 



[8] R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics 
of polymeric liquids Vol. 1 (John Wiley and Sons, New 
York, 1987). 

[9] R. G. Larson, Constitutive equations for polymer melts 
and solutions (Butterworths, Boston, 1988). 



[10] R. R. Huilgol and N. Phan-Thien, Fluid mechanics of 
viscoelasticity (Elsevier, Amsterdam, 1997). 

[11] M. Vladkov and J. L. Barrat, Macromol. Theory Simul. 
15, 252 (2006). 



